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PRIMITIVE EQUATIONS OF THE OCEAN 



E. AUDUSSE, P. DREYFUSS, B. MERLET. * 

Abstract. In this article we are interested in the derivation of efficient domain decomposition 
methods for the viscous primitive equations of the ocean. We consider the rotating 3d incompressible 
hydrostatic Navier-Stokes equations with free surface. Performing an asymptotic analysis of the 
system with respect to the Rossby number, we compute an approximated Dirichlet to Neumann 
operator and build an optimized Schwarz waveform relaxation algorithm. We establish the well- 
posedness of this algorithm and present some numerical results to illustrate the method. 
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1. Introduction. A precise knowledge of ocean parameters (velocity, tempera- 
ture...) is an essential tool to obtain climate and meteorological previsions. This task 
is nowadays of major importance and the need of global or regional simulations of the 
evolution of the ocean is strong. Moreover the large size of global simulations and the 
interaction between global and regional models require the introduction of efficient 
domain decomposition methods. 

The evolution of the ocean is commonly modelized by the use of the viscous primitive 
equations. This system is deduced from the full three dimensional incompressible 
Navier-Stokes equations with free surface with the use of the hydrostatic approxima- 
tion and of the Boussinesq hypothesis. It is implemented in all the major softwares 
that are concerned with global or/and regional simulations of ocean and/or atmo- 
sphere (we refer for example to NEMO [23j, MOM [26j or HYCOM for global models 
and ROMS [2] or MARS for regional models). The primitive equations have been stud- 
ied for twenty years and important theoretical results are now available [211 E21 H] • 
The numerical treatment of this system has been also strongly investigated [31] . But 
the key point here is to simulate global circulation on the earth for long time and/or 
with small space discretization. This type of computations can not be performed on 
a single computer in realistic CPU time and need to be parallelized. The problem is 
then to allow the different subdomains to interact in an efficient way. Another type 
of applications that are commonly investigated in the oceaonographic and/or meteo- 
rological community is to couple global and regional models in order to obtain precise 
regional previsions. The problem is also to construct an efficient interaction between 
the two models. In [3] the authors exhibit that most of the existing algorithms are not 
able to compute this kind of problem in an efficient way. We propose in this article 
to investigate these still open questions in the context of a quite recent performing 
domain decomposition method : the Schwarz waveform relaxation type algorithms. 

The development of domain decomposition techniques have known a great develop- 
ment for the last decades and our purpose is not to make an exhaustive presentation 
of these methods. We refer the reader to [27l [33] for a general presentation and we 
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restrict ourselves to the description of Schwarz waveform relaxation method. It is a 
relatively new domain decomposition technique. It has been developed for the last 
decade and has been successfully applied to different types of equations. This type 
of algorithms is the result of the interaction between classical Schwarz domain de- 
composition techniques and waveform relaxation algorithms. Its great interest is to 
be explicitly designed for evolution equations and to allow different strategies for the 
space time discretization in each subdomain. Moreover we can even consider different 
models in each subdomain without modifying the architecture of the interaction. 
The heart of the classical Schwarz method is to solve the problem on the whole domain 
thanks to an iterative procedure where a problem is solved on each subdomain by the 
use of boundary conditions that contain the information coming from the neighboring 
subdomains. It comes from the early work of Schwarz where this idea was intro- 
duced to prove the well-posedness of a Poisson problem in some nontrivial domains. 
This method is designed for stationary problems and presents two main drawbacks : 
it needs an overlapping between subdomains and it converges slowly [19j . In the last 
decade, some works have been devoted to cure these disagreements [20]. We refer to 
[To] for a complete presentation. 

The extension to time evolution problems was performed at the end of the nineties by 
Gander [51 [5] and Giladi & Keller [T3] and was denoted Schwarz waveform relaxation 
algorithms. The authors mixed the classical Schwarz approach with waveform relax- 
ation techniques developed in the context of the solutions of large system of ordinary 
differential equations [HI [17] . The exchanged quantities were of Dirichlet type. Op- 
timized Schwarz waveform relaxation methods were developed with the introduction 
of more sophisticated information to compute the interaction between the subdo- 
mains. These optimized algorithms were based on previous works [71 1151 116] about 
the derivation of absorbing boundary conditions respectively for hyperbolic, elliptic 
and incompletely parabolic equations. The same ideas were used to derive efficient 
transmission conditions between the subdomains : since the exact transparent condi- 
tions can not be implemented in general (it may lead to non-local pseudo-differential 
operators), the derivation of some approximate conditions is performed. These condi- 
tions can be optimized with respect to some free parameters which justifies the name 
of the method. The optimized Schwarz waveform relaxation method was first applied 
to the wave equation [l^ and then to the advection-diffusion equation with constant 
or variable coefficients [24] . A recent paper [11] gives the complete solution of the one 
dimensional optimization problem for constant coefficients equations. More recently 
the method has been extended to the linearized viscous shallow water equations with- 
out advection term by V. Martin [25] . Here we are interested in the application of the 
method to the system of Primitive Equations of the ocean. It leads to non-trivial new 
problems (new transmission conditions, well-posedness of the problem, convergence 
of the algorithm...) that we address in this article. 

The outline of the paper is the following : in Section [2] we write the equations and we 
precise the asymptotic regime that we consider. In Section [3] we derive an approxi- 
mated Dirichlet to Neumann operator, and define the associated Schwarz waveform 
relaxation algorithm. In Section [4] we define a weak formulation of the problem on 
the whole domain and prove that it is well-posed in the natural functional spaces. In 
Section [5] we introduce a weak formulation for the Schwarz waveform relaxation algo- 
rithm and prove that each sub-problem solved in the algorithm is well-posed. Finally 
we present some numerical results in Sections M 
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2. The set of equations. We first write the primitive equations of the ocean. 
Then we present the simphfied system from which we are able to derive efficient 
transmission conditions. 

2.1. The primitive equations of the ocean. We consider the primitive equa- 
tions of the ocean on the domain {x, y,z,t) E Rx Rx [—H{x, y), ({x, y, t)] x R+ where 
—H{x,y) denotes the topography of the ocean and ({x,y,t) denotes the altitude of 
the free surface of the ocean. The primitive equations are commonly written [5] 

dtU^ + U,, ■ V^C/^ - ,^AU,, + -f1 A U,, + -V,p = 0, (2.1) 

Po Po 

■ U,, + d,w = 0, (2.2) 

dzP = -pg, (2.3) 

p = p(z,T,5), (2.4) 

dtT + Uo-VT- i^tAT = Qt, (2.5) 

dtS + Uo-S/S~iysAS^Qs, (2.6) 

where the unknowns are the 3d- velocity ([/^, w) = (m, v, w), the pressure p, the density 
p, the temperature T and the salinity S. The parameters are the gravity g, the eddy 
viscosity v, the eddy diffusion coefficients for the tracers vt and vs and the earth 
rotation vector fi. The source terms Qt and Qs for the temperature and salinity 
model the influence of the sun, rivers and atmosphere for these tracers. 
Note that we consider here the classical but non-symmetric viscosity tensor 

dyU 
dyv 

dxW dyW 

Other form of the viscosity tensor can be found in [13j . Note also that it is possible to 
consider different viscosity coefficients in the horizontal and vertical directions [22] . 
These equations are supplemented by initial and boundary conditions. At initial time, 
we impose 

Uhi;0) = C/,,, inn, C(-,0)=C. inc., 

where the subscript letters i means "initial" . At the bottom of the ocean we impose 
a non-penetration condition and a friction law of Robin type {ah > 0) 

U^i-H)-V^{H)-w{-H) = 0, dnUti~H) + abUti-H) = 0, (2.7) 

where Ut stands for the tangential velocity and n denotes the outward normal vector 
to the bottom of the ocean. 

The free surface is transported by a kinematic boundary condition 

dtC + u^iO-v,X-MO = 0. (2.8) 

The equilibrium of the stresses at the free surface implies 

k - (P - Pa)Id] ■ ^ \ dyC \ ^ 0, (2.9) 

where Pa{x, y, t) denotes the atmospheric pressure. 
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Fig. 2.1. Schematic representation of the ocean 



2.2. A linearized hydrostatic model. In order to derive simple and efficient 
transmission conditions for the Schwarz waveform relaxation method we make some 
assumptions on this set of equations. 

First we neglect the influence of the tracers (temperature and salinity) on the density. 
Thus we suppose that the density is constant (we assume pq — 1) and we do not 
solve the equations on the tracers (|2.5p - (|2.6p . Note that these equations are classical 
advection-diffusion equations for which the optimized transmission conditions are well 
known [nilll]. 

Then we use the divergence-free condition (|2.2p and the non-penetration condition (|2.7p 
to write the vertical velocity w as a function of the horizontal velocity C/^ and we use 
the hydrostatic assumption (|2.3p to write the pressure p as a function of the water 
height C,. The remaining unknowns in the system are the horizontal velocity Uf^ and 
the water height C- The set of equations (|2.ip - (|2.3p stands 

J-H 

with 

C= J ^, and f:^2n-e,. 

The first equation is written on the initial domain R^; xR^ x [—H{x, y), ({x, y, t)]^ xR^ 
while the second one is written on R^; x Rj, x R^. We consider for simplicity a 
flat bottom and a constant atmospheric pressure. Then we linearize the problem 
around a constant state which corresponds to a horizontal velocity Uq = (Mo,fo) and 
a horizontal free surface located at z = 0. It follows that the water height C is a small 
perturbation. In the sequel denotes the perturbation on the horizontal velocity. 
The linearized problem stands 



dtU^ + U,-W„U^~ J. AC/, + fCU„ + .gV^C = 0, 
dtC + HV, • [/, + U„ ■ V,C - 0, 



(2.10) 
(2.11) 
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where 



u 



= I :=- U^dz, 



V H 



H 



denotes the mean horizontal velocity of the flow. This mean velocity is called barotropic 
velocity by the oceanographic community while the deviation Uj^ — Uf^ is called baro- 
clinic velocity [5]. 

Note that the first equation (|2.f Op is now written in the fixed domain R^^^ x Rj, x 
[-77,0], xR+. 

The associated boundary conditions are 

dMz^O) = 0, {d,U^ + abUh}{z^-H) = 0, (2.12) 

where the boundary condition at z = is deduced from the equilibrium of the stresses 
at the free surface (|2.9p . Indeed, with the help of the linearization procedure, we first 
deduce that at first order dzU,^{z = () ~0 and then dzUf^{z = 0) ~ since we assume 
that C, is small. 

In order to derive the transmission conditions we assume = in the sequel. How- 
ever in the definition of the Schwarz waveform relaxation algorithm and in the nu- 
merical simulations, the condition «{, > will be supported. 

2.3. Dimensionless system. We choose characteristic horizontal and vertical 
lengths (denoted L and H respectively) and velocity U of the problem. We introduce 
the dimensionless quantities 

{x,y) = L{i,y), t={L/U)i, 
C = HC, z = Hi, 

[/, = Ul\, u, - uu,. 

The spatial domains of computation are Q — R^. x R^ x (— 1, 0), for the momentum 
equation and uj — Tlx x Ry for the continuity equation. We study both equations in 
the time interval [0,T], where T > is fixed. 

Dropping the "~" for a better readability, the system in dimensionless variables stands 

dtU,, + Uo ■ V,C/, - ;^A,C/, - -^d^U^ + -CU^ + -^V;,C = 0, (2.13) 

d,U^ix,y,0,t) = d,Uf^ix,y,^l,t) ^ 0, (2.14) 

[/,.(.,0) = (2.15) 
5tC + f/o-V,C + V,-C7;, = 0, (2.16) 
C(-,0) = 0- (2.17) 

We have introduced the characteristic quantities, 

- e = U/{fL) the Rossby number, 

- Re = UL/v the horizontal Reynolds number, 

- Re' = H^/L^Re the vertical Reynolds number, 

- Fr = U/y/gH the Froude number. 



We choose to exhibit the Rossby number as a small parameter since we are interested 
in long-time oceanographic circulation for which the Rossby number is typically of 
magnitude 10^^. The values of Reynolds and Froude numbers vary with respect to 
the turbulent processes and to the depth of the area that is considered respectively. 
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3. The optimized Schwarz waveform relaxation algorithm. We are now 

interested in finding efficient transmission conditions for equations (|2.13p - (|2.17p . We 
first present tlie Scliwarz waveform relaxation method. Tiien we derive tlie relations 
satisfied by the optimal transmission conditions. Since we are not able to solve analyt- 
ically these equations, we perform an asymptotic analysis with respect to the Rossby 
number e in order to derive some approximated transmission conditions. Finally we 
present the related optimized Schwarz waveform relaxation algorithm. 

3.1. The Schwarz waveform relaxation method. The heart of the method 
is the following. We first divide the computational domain into an arbitrary number 
of subdomains. Then we solve each sub-problem independently for the whole time 
interval. The interactions between neighboring subdomains are entirely contained in 
the boundary conditions. An iterative procedure is considered until a prescribed preci- 
sion is reached. The advantages of the method are clear : the parallelization is almost 
optimal: at each step the sub-problems are solved independently, so the space-time 
discretization strategies (or even the models...) can be chosen independently on each 
subdomain. Moreover at the end of each step only a small amount of informations are 
exchanged. The main drawback is related to the needed number of iterations : the 
method is efficient if it converges quickly (in two or three iterations typically). This 
requirement needs the derivation of efficient transmission conditions. 
In the sequel we consider for simplicity two subdomains but the method extends to 
an arbitrary number of subdomains. 

We begin with some notations. First we introduce the left and right spatial sub- 
domains Cl^ and defined by: 

n- :^ (-C3o,0), X Rj, X (-1,0),, n+ := (0, -l-w), x Rj, x (-1, 0)„ 

and their interface 

r = {0}, X Rj, X (-1,0). ~ Ry X (-1,0).. 

We also introduce the domains uj^ := ±(0, -l-oo)a; x Ry for the unknowns that do not 
depend on the z variable, and their interface 7 :— {0}x x Ry — Ry 
Let D be some spatial open domain and T > be a given real number. Then we will 
write Dx to denote the cylindrical domain Dx := D x (0,T). 

We denote by PE the set of equations ([^T^ . ([^1^ and (PTB)) and X := (C/„,C) 
stands for the solution of this system with associated initial data Xi := (C/^ ^, Q). 
Then the Schwarz waveform relaxation algorithm is defined as follows: 



' PE(X!!+i) ^ 





on 17^, 




' PE(X!f+i) = 





on 


< X!!+i(-,0) = 


x^ 


on 0", 


< 


Xl+\-,Q) = 


x^ 


on fl+ 




B-Xl 


on Ft, 




. B+XI+' = 


B+X": 


on Ft 



where the operators B± contain the transmission conditions. 

In the classical Schwarz waveform relaxation algorithm [8l [6] , the transmitted quan- 
tities are of Dirichlet type and the operators B± are thus chosen to be the identity 
operator. Note that in this case an overlap is needed in the definition of the subdo- 
mains. 
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In the sequel we are interested in deriving more efficient transmission conditions. In 
order to reach such a goal we will first describe the general method to obtain optimal 
transmission conditions. The transmission conditions are said to be optimal if the 
algorithm converges in two iterations to the solution of the initial problem. These op- 
timal transmission conditions involve the Dirichlet to Neumann operator associated 
to PE on the subdomains fl^. Here we will see that we are not able to obtain an 
explicit formulation for these optimal conditions. Anyway these optimal boundary 
conditions are not local and consequently too expensive to be useful from a numerical 
point of view. 

Recent methods have been developed recently in order to approximate these optimal 
conditions by analytical or numerical means — see the review paper [TU] for elliptic 
problems and for parabolic evolution equations. 

Here we will perform an asymptotic analysis of the system with respect to the Rossby 
number e in order to deduce a set of approximated and efficient transmission condi- 
tions. This strategy has been initiated in |25| for the shallow water equation without 
advection term. In our case, it turns out that these approximate transmission condi- 
tions lie in a two parameter family of boundary conditions. In Section [6] we optimize 
numerically the transmission conditions in this two parameter family. 

Let us first describe in a formal setting the ideal case of optimal transmission condi- 
tions for the Schwarz waveform relaxation algorithm (j3.ip . 

We consider the case uq > 0. The case mq < is deduced by applying the symmetry 
"a;' = — x" . Integrating the linearized Primitive equations PE on a subdomain, we 
see that the flux of the unknown ([/^, through the interface F is given by 

■^d^U^ - ^oUh - ( ) ' + " 

Using this flux as a Neumann operator, we define the Dirichlet to Neumann operators 
as follows. Consider a Dirichlet data Xi, = (?7^ j,, Cb), we set 

DN'^^X, := - uoU, - ( ^ c 

where X = {Uf^X) solves 

PE{X) =0 on n+ 
X(-,0)= onn+, 
X = Xb on Tt- 

Symmetrically, consider a Dirichlet data C/^ 5, we set 



\ DNiXt J 
where X — {Uf^X) solves 




PE(X) =0 on r?^, 
X{-,0)= onn-, 
{/, = C7, on Ft. 
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Notice that since we consider the case uq > 0, the continuity equation (|2.16p the 
boundary condition is relevant only in the subdomain 57 J. This is why in the later 
case we do not have to prescribe a boundary condition for C. 



Once these Dirichlet to Neumann operators are defined we can introduce the optimal 
transmission conditions 



B+X = 



uoC + u- DNiX 



DN^'^X 



dnI'^x \ 



I 



(3.2) 
(3.3) 



Proposition 3.1. With this particular choice of transmission operators B±, the 
algorithm (j3.ip converges in two iterations. 



Proof. By linearity, we may assume that the exact solution is {Xi = 0). At the 
initial step the solutions on each subdomain do not satisfy any particular property. 
But the first iterate solves the primitive equations with vanishing initial data. It 
follows from the very definition of the operators DN± that in the definition of the 
second iterate, the right hand sides of the transmission conditions vanish for both 
sub-problems. We deduce that this second iterate vanish: the algorithm converges in 
two steps. □ 

The operators (|3.2p(|3.3p being non-local pseudo-differential operator, they are not well 
suited for numerical implementation. Our strategy is to approximate these operators 
by numerically cheap operators. Of course the two-step convergence property will be 
lost. The quality of the approxiamation will be measured through the convergence 
rate of the algorithm. From the structure of (|3.2|)()3.3|) . we choose to write B± as 
perturbations of the natural operators transmitted through the interface: 



B-X = 



B+X 



V 



mqC + u — S'lX 



(3.4) 
(3.5) 



where S^'' and are pseudo-differential operators that will approximate the Dirich- 
let to Neumann operators. 

Let us finally remark that the differences in the expression of the two transmission 
operators B± are due to the sign uq > 0. Since B- contains the information that 
is transmitted from fl^ to fl^ it is constructed on three boundary values (velocities 
and water height) but it has to transmit only two boundary conditions for momen- 
tum equations (|2.13p . On the contrary B^ is constructed on two boundary values 
(velocities) but has to send three boundary conditions (for momentum and continuity 
equations). 

In the next subsections we will identify optimal and approximated transmission op- 
erators. To carry out the computation of the Dirichlet to Neumann operators we 
perform Fourier-Laplace transforms. 
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3.2. Laplace-Fourier transform of the primitive equations. We perform 
on the set of primitive equations (|2.13p - (|2.16p a Fourier transform in the y variable 
and a Laplace transform in time. The dual variables are respectively denoted 77 G R 
and s = a + ir Cz C. The real part a is assumed to be strictly positive. We obtain in 
each subdomain the same set of differential equations 

{s + uodx + ir]Vo} ( + dxU + irjv = 0. 



In the z direction we introduce the eigenmodes of the operator —d]^ on (—1,0) with 
homogeneous Neumann boundary conditions (|2.14p 

e„(z) :— Q;„cos(/i„z) with :— nir; ao :— 1 and a„ := V2 if n > 0. 

Then we search for the solution on the form 

00 

;7Ua;,^) = 5]?/r(a^)e„(^). 

Note that we obviously obtain Uf^ = C/°. It means that the first vertical mode C/° rep- 
resents the barotropic velocity while the sum of the other ones denotes the baroclinic 
deviation. 

The barotropic mode is coupled with the water height and it is the solution of the 
following system of three ordinary differential equations, 

- ^e''^'^ + + + ^^^^ + + ^° + ( ) ^ ^'-'^ 

uod^C +{s + irjvo) C + 8.^11^ + ir/v^ = 0.(3.7) 

This last system is exactly the Laplace-Fourier transform of the so-called linearized 
viscous shallow water equations [28]. 

For the other vertical modes we have a set of two coupled reaction advection dif- 
fusion equations, 

- ^9^UJ: + uod^UJi -I- |s + iiivo + ^rf + + Jsj UJ^ - 0. (3.8) 

3.3. Optimal transmission conditions for the baroclinic modes. The 

derivation of optimal transmission conditions for an advection diffusion equation was 
performed in [21] . Here we are interested in the set of coupled reaction advection dif- 
fusion equations (|3.8p . The baroclinic modes are not coupled with the evolution of the 
water height. Hence for these modes the transmission operators have two components 
and will be searched on the form 

Bl = ( T^^d,,U^±uoU„-S^- Y (3.9) 
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We search for the solution of system p.Sp as a sum of exponentials x i-^ e^^ . Plugging 
this ansatz in the system, we obtain that e^^ solves (|3.8p if and only if A is a root of 
the determinant of the matrix M„(A) := 

This determinant is a polynomial of degree four in A and we can compute its four 
roots 

Ai^+ ^ {no + ^) , Xr := ^ (up - , (3.10) 

where 

AS ■-ul + ±i^f- + § + s + ^r,.,±l). (3.11) 

Every solution A!['^ is associated with a one dimensional kernel generated by the 
vector 0!^'^ defined by 

Since the solutions must vanish at infinity we search for solutions in Q~ on the form 

Ul_{x) = a'l:+e^+^''(f>'^+ +al-+e^-^ ^ = • exp (xA"'+) • a"'+, (3.12) 

where 



:=(_',;). A- := ( f ) . := ( ) . (3.13) 

In f2+ we search for the solution on the form 

Ul+{x) = a"'"e^+"^0"'" +a"'"e^""^(/)"'" = • exp (xA"'-) • a"'". (3.14) 

It follows from relations ()3.12|) and (I3.14|) that 

dxUh.^{x) ^ • A"'=^ • exp (xA"'±) • = • A"'± • • ^7^,^. (3.15) 

We can now define the operator 5^'" in (j3.9p in order to derive an optimal algorithm. 
This is done through its Laplace-Fourier symbol: 

5^-" =F-i-$"-±A"'± ±uo Id. (3.16) 

3.4. Approximate transmission conditions for baroclinic modes. Since 
we want to construct an efficient but simple Schwarz waveform relaxation algorithm 
we will derive approximated transmission conditions by considering an asymptotic 
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analysis of the results of the previous subsection. 

The definition (|3.1ip of A!^ leads to the expansion iSJ — ^^^'"app + 0{y/e) with 




(3.17) 



Note that the approximated operator ()3.17p does not depend on n. Consequently 
the related approximated transmission operators l|3.9p can be applied to the whole 
baroclinic velocity, i.e. to the sum of the baroclinic modes. 

3.5. Approximate transmission conditions for the barotropic mode. 

The derivation of optimal transmission conditions for the linearized viscous shallow 
water equations without advection term was performed in 25J. Here we are inter- 
ested in the linearized viscous shallow water equations (I3.6|) - p.7|) . The transmission 
operators will be searched on the form ()3.4|) - (|3.5|) . 



As for the baroclinic modes we search for the solution of system (|3.6p - (|3.7l) as a 
sum of exponentials e^^. Here A has to be a root of the determinant of the matrix 
Mo(A) defined by 

^-^ + z.oA + . + £ + ^,.o -J X/Fr^ ^ 

1 A^ T]'^ i-q 

\ A if] s + uqX + i?]VQ j 

This determinant is a polynomial of degree five which does not admit a trivial de- 
composition. Hence it is not possible to derive an explicit formula for the solutions 
of p.6|) - ()3.7p . Consequently, we are not able to obtain an explicit form for the optimal 
transmission conditions for the barotropic mode, even in Fourier-Laplace variables. 
In order to derive approximated transmission conditions we use the fact that the 
Rossby number is a small parameter to compute approximated values of the roots of 
the determinant of Mo(A). The related approximated transmission conditions will be 
coherent with the results of the previous subsection for the baroclinic modes. 
Since uq is positive we first notice that three roots (|3.18p - (|3.19p have a negative real 
part and two roots p.20p have a positive real part. The negative roots will be denoted 
A^" and Ag. The positive ones will be denoted A^^. The notations for the related 
quantities that we introduce later are coherent with the previous ones (|3.13p . As 
above, we search for the solution in 17^ on the form 

X_{x) = a'^+e^+^^0'}:+ -I- a'^+e^-^-0!^+ • exp {xA'>'+) ■ a°-+. 

In f7+ we search for the solution on the form 



=: • exp (xA°'-) • 
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We compute the following approximations for the roots of the determinant of Mq{X): 

s + ir]Vo 



Aq — — 



Uo 



g _ \/±iRe I Re uq 
X± = ^ h 



Re \ 



2 AFr'^uo j 



Q _|_ V ±iRe ( Re uq 



Re 



2 AFr'^uo I 



0{V~e). 



(3.18) 
(3.19) 

(3.20) 



The associated kernel is always one dimensional and spanned by: 



$0 





I +0{6% 



iV2R( 
AFr"^ 



$1- = 



4 uoFr^ 2 



Uo - 



4Fr2 
iV2Re 



\/2 V2 1 - i 



:((±Uo + Wo)?? + s)| 



0{e) 



V 



4 Moi^r2 2 Vife 



:((-uo + TOo)'? + s) > 



0{e). 



I 



As in the baroclinic modes case, we compute the approximated transmission operators 
in Laplace-Fourier variables by 



2,3 



Mo 

Uo 



l-<£> L -I 



aPP~ Re 

where Mi i, denotes the first 2x3 matrix extracted from the 3x3 matrix M. It leads 
to the following Laplace-Fourier symbols 
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and 



^0 

'^+, app 



/ V2 



'Ree 
V2 



Mo 



1 

Fr'^UQ 
1 



n/2 1 

/lie 2Fr2uo 

v/2 



V 











0(\/i).(3.22) 



By using relations p.l7p . (|3.2ip and p.22p . we notice that 



1 



1 



S. 



±, app 



2,2 



Au.n , 

'^±, app + 



2Fr^UQ 
1 







It follows that a part of the transmission conditions will be applied to the whole 
velocity (sum of baroclinic and barotropic modes) while a second part will be applied 
only to the barotropic mode. The first part corresponds to the operator 5± "app- The 

second one corresponds to the remaining terms in the operator 5^''^app- 

3.6. The optimized Schwarz waveform relaxation algorithm. Thanks to 
the computed approximated operators ()3.17p . (|3.21|) and ()3.22p we can now derive 
an approximated Schwarz waveform relaxation algorithm for the linearized primitive 
equations (|TT3| - ((2l7|) . 



Since the computed operators p.l7p . p.2ip and (|3.22p do not depend neither on 
the Fourier variable 77 nor on the Laplace variable s the related operators in the real 
space are identical to their Laplace- Fourier symbols. It follows that the approximated 
transmission operators B± p.4p - (|3.5p have the following form 



/ 1 



Re 



■dxU 



Re 



B+X = 



1 

Re 



V2 


— u — 


V2v 


2^ Ree 


2^/Ree 


V2 




V2u 


2y/Ree 




2VRee 


C , / 


V2 





:-v/2 \ 



2Fr^uo 



(3.23) 



+ v/2 \ 




AFr^uo 



V 



,(3.24) 



for which we recall that u and v represent the mean-values with respect to the z 
variable of the velocities u and v. 

Note that by replacing the first component {B^X)i by the linear combination 



(B+X), - l/{Fr^u^){B+X\ 
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we replace p.24p by the equivalent transmission conditions 



/ 1 



BZX = 



Re 



■dxV 



v/2 \ 



\/2u u 
2^/Rle^ 4:Fr^uo 



V 



UqC + U 



(3.25) 



In the sequel we use p.25p rather than (I3.24p and we drop the superscripts "~" . 
Next, we remark that the transmission conditions p.23pp.25p are a particular case 
of the generalized transmission conditions 



1 



Re 



-dM, 



Uo 



where 



B+X = 



A := 



UqC + U 



1 -1 
1 1 



B := 



1 -1/2 
-1/2 



The original transmission operators p.23p - p.25p correspond to the choice 

1 1 



'2Re' 



/3 = 



(3.26) 
(3.27) 

(3.28) 
(3.29) 



2Fr'^UQ 

Notice that B-X and (S+X)(i 2) do not depend on the water height so we may 
rewrite B-X = B^J^U,^ and B+X = \B+'^U^,B%X) as 



Let us emphasize the identity: 



bI-u^ + b'!-u,,^2^au,,. 



B^X := uoC + u. (3.30) 



(3.31) 



This relation will be useful both for defining a weak formulation of the algorithm in 
Section [5] and for the numerical implementation of this algorithm in Section \6\ 
Finally the Schwarz waveform relaxation algorithm p.l|) writes 



f PE(X!!+i) 











' PE{Xl+^) = 





on 







on 




















Xl+\;0) = 


x^ 


on 




X, 


on 




< 


















bI'^u^X'^ 


bI-ui_ 


on 


Ft, 


on 


Ft, 






B'lX'l 














on 


7t- 



(3.32) 
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where the operators B^'' , are defined by equahties (|3.28|) (|3.30p and where a and 
f3 are free parameters. 

These generahzed transmission conditions can now be optimized with respect to the 
two parameters a and f3. In the case of a one dimensional reaction advection diffusion 
equation this optimization problem has been solved analytically (see [H]). Here, we 
will present a numerical procedure in Section [S] 

4. Well-posedness of the linearized Primitive Equations. In the previous 
sections we have performed formal computations on the linearized Primitive Equa- 
tions leading to the construction of the Schwarz waveform relaxation algorithm (|3.32p . 
The aim of this section is to be more precise: we will define a weak formulation of 
the system (|2.13p - (|2.17p and then prove that this system is well-posed in the natural 
spaces associated to this weak formulation. 



From now on we relax the boundary condition on the bottom, i.e. we assume af, > 
instead of at = 0. Moreover, in order to prepare the study of the well posedness 
of the algorithm (|3.ip in the next section, we consider non-homogeneous right-hand 
sides Y = {Fi,F2,f) = Y{x,y, z,t). The system of linearized primitive equations 
PE(X) = Y writes 

1 1 , 1 1 1 





F 


in ^It, 


(4.1) 


dzU^{x,y,0,t) = 





on LOT, 


(4.2) 


-9zUf,{x, y, -1, t) + abU,,{x, y, -1, t) = 





on ujt, 


(4.3) 




/ 


in lut- 


(4.4) 



We supplement this system with the initial conditions 

[/,(•, 0)= inn, (4.5) 

C(-,0) = inw. (4.6) 

Note that if we consider that the water height is given, the system (|4.ip - (|4.3p . (|4.5p 
with unknown [/^ is a classical linear parabolic problem. On the other hand if we 
consider that the mean horizontal velocities JJ/j are given then C solves the linear 
transport problem with source term ()4.4p . (|4.6p . 

We will proceed as follows: first we recall the classical weak formulations both for 
the parabolic problem (with prescribed water height) and for the transport equations 
(with prescribed velocity). These two problems define two maps Si : C, ^ Uy^ and 
S2 '■ 1-^ (. Finally we define the weak solutions of the Primitive Equations to be 
the fixed points of the map r : (J7^,C) '-^ i^iiO, ^2{Uf^)) and conclude by proving 
the existence of a unique fixed point. 

Let us first introduce some functional spaces and some notations. We will work 
with initial data and right hand sides satisfying 

if :=L2(f],R2), a^L^Lo), 

Fe L2(0,r;V'), f eL\0,T;LHLu)), 

where V' is the topological dual of V := H^{fl, R^). 
The weak solutions will satisfy 

{7„e C{[0,T],H)n L\0,T;V), ( £ C ([0, T], L2(c^)) n C7(R,, (R^ x (0, T))) . 
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We will need the following bilinear forms: 

<u,v) := i- (v,[7,v,y)^^ + id.u,d.y),,^ + ^iu,v)^_,,^ 

+ ^iCU,V)n^ + iUo-VU,V)n^, (4.7) 
ciCV):=^{Ce.,d.,V)^^. (4.8) 

where w_i :— Hx x Rj, x { — 1}^ and {U, V)j: denotes the scalar product on S. 

Assuming that we have a strong solution, and taking the scalar product of equa- 
tion (|4.ip with V e Vlfl X (0,T),R^), we obtain (after integrating by parts) the 
following weak formulation for the equations governing the horizontal velocities: 

Vy el?(f7x (0,T),R2), {^tU^,,V)^^+a(U„,V) = c{C,V) + {F,V). (4.9) 

We now state 

Definition 4.1. Let F e L'^{0,T;V') and ( e L'^{lot), we say that U,^ e 
L^{0,T;V) is a weak solution of the system (j4.ip if (|4.9p holds. 

Proposition 4.2. Let C/^_, e L^{n), F e L^{0,T;V') and ( e L^{ujt), there 
exists a unique weak solution Uf^ € C{[0,T]; H) C] L^{0,T;V) of (j4.ip satisfying the 
initial condition (j4.5|) . Moreover, we have the energy inequality 




< l\\UhAl + j\{F,U^){s) + idxuXUs)}ds. (4.10) 



Proof. The method is classical and we only sketch the proof. We obtain the 
existence of a solution satisfying (|4.10p by the Galerkin method. Let (Em) be an 
increasing sequence of finite dimensional sub-spaces of V such that Ui?m is dense in 
V. For every m, there exists G C°°([0, T], i^^) such that (|4.9p holds for every 
V e I'(0, T; En) - we only have to solve a finite system of linear ordinary differential 
equations. 

Using Um X l[o,t] as a test function in (|4.9p . we conclude that fJ™ satisfies (I4.10p . 
So using the Cauchy Schwarz inequality and the Gronwall Lemma, we see that the 
sequence {Um) is uniformly bounded in i^(0, T; V). 

Now from the weak formulation, we deduce that [dtUm) is bounded in L'^{0,T;V'), 
thus by Aubin-Lions Lemma, [Um) is compact in C{[0,T],H). Extracting a subse- 
quence we obtain a solution satisfying (|4.10p . 

Regularizing in time and using the weak formulation, we see that any solution satis- 
fies (|4.10p and uniqueness follows by the energy method. □ 

Let us turn our attention to the equations ()4.4p . ()4.6p governing the evolution of the 
water height ^. This is a linear transport equation with constant coefficients and a 
source term. Assuming that is a strong solution, multiplying (j4.4p by a test function 
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X G "Diij X [0, T)), integrating on wy, integrating by parts in space and time and then 
using the initial condition (|4.6p . we obtain 

Vxel?(^x [0,T)), 

-(C,{at + c/o.v4xU- (C.,x(-,o))^ + (/-v,.I7,„xU. (4.11) 

Definition 4.3. Let f € L^itot), C/^ G L2(0,T; V) a^d Q e L'^{lo), we say that 
C e L'^{ll!t) is a weak solution of the system (|4.4p . Ij4.6p «/ (|4.1ip holds. 

Remark that the test function does not necessarily vanish at time and that the 
initial data is prescribed by the weak formulation. 

Proposition 4.4. Let f e L^{ujt), U^^ e L^{0,T;V) and Q e L^{Vl). There 
exists a unique weak solution C, G L'^{ujt) of (|4.4p . (j4.6p . Moreover this solution is 
given by the characteristic formula: 

C{x,y,t) ^ dix - uot,y - vot) + if -Vh ■Uh){x - uos,y - vos,t- s)ds. (4.12) 

Jo 

This solution lies in C ([0, T]; _L^(w)) n C (Rx] L'^{'R.y x (OjT))) and satisfies the fol- 
lowing estimates for every t G [0, T] and every a; G R, 

||C(-,-,OIU < / ||/-V,-C7j|^(s)ds, (4.13) 
Jo 

\\a^,;-)K<-(\m^+ [ \\f-Vh-U,,Us)ds). (4.14) 
"^0 \ Jo J 

Proof. First, notice that the estimates (|4.13p (|4.14p are direct consequences of the 
characteristic formula (|4.12p . 

Next, remark that if the data V/jt/^j, / and (i are sufficiently smooth then the function 
C given by the formula (|4.12p solves (|4.1ip . Hence we obtain the existence of a solution 
of (|4.1ip by density. 

For the uniqueness, by linearity we may assume that the data V/iC//j, / and Q vanish. 
Then let V' S 2?(a;) and p G T>{[0,T)) and define the test function x by xi^iVii) •= 
—'ip{x — uot, y — VQt) p{s)ds, so that: 

9tX + Uo-Wx = - Uot, y - VQt)p{t), 

and (|4?TT|) yields 

= / C{x,t)p{t)ilj{x - UQt,y - vot) = / Cix + UQt,y + vot,t)p{t)ilj{x,y). 

Since this is true for every (^/;, p) G T^iyj) x I?([0, T)), we have C = on lot- Q 

Finally, we define the notion of weak solution for the linearized primitive equations. 

Definition 4.5. Let Y = {F, f) G L^{nT,V') x L^{ujt) and Xi = (U^ .Xi) & 
L'^iQ) X L'^{uj). We say that X = (J7^,C) ^ C{Q,T;H) x L'^{ujt) is a weak solution 
of (HTT]) - ((iTB)) if the weak formulations (j?^ and pTTTjl hold and ifUf^{-,0) = L/^^. 
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Theorem 4.1. Let Y = {FJ) e L^{nT,V') x L^{lut) and X, = C») G 

L^(f2) X L'^ibj). There exists a unique weak solution X = {Ui^,C) G (C(0,T;_ff) n 
L^O,T;V))x L^{ujt) of l^-l^. 

Proof. The right hand side Y and the initial data Xi being fixed, Proposition [121 
and Proposition 14.41 define two maps 

Si: L\coT)^CiO,T;H)nL^O,T;V), ( ^ U^, 

and 

^2: L2(0,r,V)-.C([0,r];i2(,^))nC(R,;L2(Rx (o,r))), ^ c 

Denoting by T the affine mapping (f//j,C) {Si{C), S2{Uf^j), the application .'C is a 
weak solution of (I4.ip - (|4.6p if and only if it is a fixed point of T in 

£r^ := L^{O.T;V) x C{[Q,T], L^{uj)). 

Let Xi,X2 e £t and let {U^,0 ■= - ^2 and (;7;,,C) := T{Xi) - r(X2), by 
linearity, using (|4.10p . we get for < i < T, 




By Young inequality, we may absorb the term in VUf^ in the left hand side and get: 




(4.15) 



for some k > 0. Now (|4.13p and the Cauchy Schwarz inequality yield 

WCWlit) < t f \\VUjUs)ds, forO<t<r. (4.16) 
Jo 

Finally, inequalities (|4.15p (|4.16p imply that, for T' G (0, T] small enough, the mapping 
T is strictly contracting in £t' yielding the existence of a unique fixed point of T in 
£t'- Repeating the argument on the intervals [T\2T'], [2T',3T'], ... we obtain the 
result on [0,r]. □ 

5. Weak formulation and well-posedness of the Schwarz waveform re- 
laxation algorithm. We study in this section the well-posedness of the algorithm ()3.32p . 
First, we will define weak formulations for the two sub-problems and prove that they 
are well-posed. We will pay a particular attention to the weak form of the transmis- 
sion conditions. In particular we will establish that the solutions XJ"*"^ of the n*^ 
step of the algorithm (|3.32p are in the right spaces, allowing the construction of the 
transmission conditions for the next step. 

As in the previous section, we also consider non-homogeneous right-hand sides Y = 
{F,f). Every step of the algorithm may be split in the two following sub-problems. 
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First in the domain {x < 0}, we search for a solution :— X_ — (C/,j ,C_) 

solving the initial and boundary value parabolic problem, 

+ ^0- V. - i^A, - + icj U,^_ + ^V,C- = F in , 

-d:,Uf^_{x,y,-l,t) + ahUf^_{x,y,-l,t) = 0, 

dzUf^^{x,y,0,t) = onw^, 

i3^''f/^ _ = ^^-UZ+ on Ft, 



(5.1) 



(5.2) 



and the transport problem, 

C+'(-,o)= 

In the right subdomain {x > 0} we search for a solution := X+ = {U/-^ +iC+) 

solving the initial and boundary value parabolic problem, 

dzU,^_^_{x,y,0,t) = onw+, 
S^-- = Bl-Ul_ on Ft, 

f^/..+ (-,0) = U^^, inn+, 

(5.3) 

and the transport problem with entering characteristics on the boundary jt, 

C+(-,0) = in^+. 

To prove that these two sub-problems are well-posed, we proceed as in Section [H 
First we study the parabolic problems with prescribed water heights: we introduce 
a weak formulation for these problems and prove that they are well-posed. Then we 
study the transport equations, introduce their weak formulations and establish their 
well-posedness. Finally, the solutions of the coupled parabolic-transport problems are 
obtained via a fixed point method. 

As in SectionlH the initial data Xi{Uf^ ^,Ci) satisfy C/^ ^ & H, Q & ^^(w). We choose 
right hand sides Y = {F, f) in L^{0,T; H) x L'^{ujt). (In section^ we only assumed 
F e L^(0,T; V), but here this choice would cause difficulties at the interface). We 
will search for weak solutions X± — {Uf^ ^,(±) in the spaces, 

C/;,,±eC([0,r],H±) n L^iO,T;V^), (5.5) 

C±eC([o,r],L2(c^±)) n c(r±,„l2 (r^ x (o,r)t)) , (5.6) 
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withi7± ;= L2(r2±,R2) and V± i/i(17±,R2). 

5.1. The parabolic problems. Let us define the weak-formulation for tlie 
parabolic problems (jS.ip and (j5.3p . First we introduce the bilinear forms and 

+ i(CC/,l/)^± + ({/o-VC/,l/)^±, (5.7) 
^^(C-^) = ^ {Ce.,d.V)^^ ± , (5.8) 

where lu'^i :— x Ry x {— l}z- 

Next, taking the scalar product of the first equation of (|5.1|) or (|5.3p with some test 
map y e ^(ITi X (0,T),R2), we obtain: 

(9tf/,,±, y),,±,(o,T) + ci^iUH,±,V) = c^{C±,V) T ]^(^-C^^,±' ^)r + ^)o±- 
Then, using the transmission conditions to express dxU/^ ^ on F, we get 

(9tC/^,±,y)^±,(o,^) +a±(C/^,±,y) + fe±(C/,^±,^) 

= c^{c±,v) + {bI- uj:^^,v)^ + {f, y )^± . 

with 

5±(C/,y) := ±^(C/,^)r + ^(AC/,y)rT/3(SC7,y)r. (5.9) 
2 ve 



We are still not satisfied with this weak formulation. Indeed, the knowledge of dxUJ^^ 
on the boundary F x (0, T) is needed for defining the term {B^'^U]^ ^, V)r in the right 
hand side. Unfortunately, (|5.5p only gives: d^UJ^-^ G ^^(il^ x (0,r)) which is not 
sufficient to define a trace. To overcome this difficulty, we use relation p.3ip to define 
recursively the terms (B^'^UJ^ ^,V)r- Indeed, for strong solutions, we have on F^ 

B-^U,^, ^ -Bl^U,^, + 2^AU,^, ^ -Bl^Ul, + 2^AU,^,. 

Thus, identifying S±'^t/^_^ with a distribution B^ G L'^{0,T;W), where W denotes 
the space _ff^/2(F, R^) ; we obtain a weak formulation of the algorithm for the hori- 
zontal velocities: 

Definition 5.1. Assuming that the functions C± = C±^"'^ ^'"c known, the weak 
formulation of the parabolic part (|5.ip and (|5.3p of the algorithm (j3.32p are defined 
as follows: 

For the first step, we choose 

Ble L^iO,T;W') (5.10) 
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Then for n > 0, the horizontal velocity is defined by UJ^^ — Uf^ ^ where U/^ solves 

= c±(C±, V) + {Bl , V)^^ + (F, (5.11) 

where , , and 6^ are defined in (j5.7p — (j5.9p . Once UJ^^ is known, we can define 
the boundary conditions for the next step in the opposite domain by 

-Bl + 2^AU;:X\r.- (5-12) 

Notice that assuming that the maps f/^^^ satisfy (|5.5p then their traces on Tt are 
weh defined in L^(0, T; W) C i^(0, T; W). Consequently, the transmission conditions 
defined recursively by ((57n)l stay in the space 2.^(0, T; W). 

Proposition 5.2. Lei t/^^^^ e H, F e L^{0, T- H), B''l e L^{Q, T; W'±) and C± (= 
C±^^) satisfying (|5.6p . T/ien i/iere exists a unique UJ^^ = twii/i regularity (j5.5p 

satisfying (|5.1ip and t/ie initial condition Ul^^{0) = Uj^^ on Q,±. Moreover, we have 
the energy inequality 




Proof. We proceed as in the proof of Proposition 14.21 we apply the Galerkin 
method. Here we only check that the a priori inequality (|5.13p is sufficient for applying 
this method. In order to bound the quadratic terms in the left hand side of (|5.13p 
and the last term in the right hand side, we will use the inequality 

\\U\\l < 2\\U\\n±\\d,U\\n±, 

valid for U E V±. (To prove it, write |C/(0, y, = 2 j"^{d^U ■ U){x', y, z) dx' inte- 
grate on Rj^ X (—1, 0)2 and use the Cauchy-Schwarz inequality). From this inequality, 
the Cauchy-Schwarz inequality, the Young inequality and the fact that the trace on 
r defines a continuous embedding 11 : — > W, we see that (|5.13p implies 




for some k > 0. Taking a Galerkin sequence {Um) associated to (|5.11l) . the elements 
of this sequence satisfy (j5.13p and then inequality (|5.14p and the Gronwall Lemma 
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imply that this sequence is bounded in L^{0,T;V). Extracting a subsequence (as in 
Proposition 14. 2p we obtain a solution of (|5.1ip . 

Then using the weak formulation satisfied by Um we see that (dtUm) is bounded 
in L^(0,T;V") and from Aubin-Lions Lemma (see e.g. [30j), the sequence (Um) is 
compact i^(0, T; iJ*(il^, R^)) for s < 1. Thus we may let m tend to oo in the 
quadratic boundary terms in the left hand side of (|5.13p . 
The uniqueness follows from (|5.14|) and Gronwall Lemma. □ 

5.2. The transport equations. We now consider that the velocities — 
UJ^'^ are known and study the transport problems (|5.2p and (|5.4[) . We begin with the 
domain {a; < 0}. Proceeding exactly as in Section 31 we obtain that a strong solution 
of Problem (|5.2p satisfies 

VxGi?(c^- X [o,r)), 

-(C-,{a, + C/o-V4xL- = {Q,x{;0)h-+if~^h-Uh.-,xh-- (5.15) 

Definition 5.3. Let f e L^iuJt), U,^ _{^ U'^+_}) € L^{0,T;V-) and Q e L^{lu). 
We say that ^_ £ L'^(uj^) is a weak solution of Problem (j5.2p if (j5.15p holds. 

The following result is proved exactly as Proposition 14.41 

Proposition 5.4. Let f e L'^icut), C/,"t^ e L^{0,T;V-) and Q e L^{uj). There 

exists a unique weak solution C"^^ = C- ^ L'^{uj^) of ()5.2p . Moreover this solution 
is explicitly given by the formula: 

C-{x,y,t) = Q{x-uot,y-vot)+ / (f -S/ h- Uh^_){x - uos,y - vos,t - s)ds. (5.16) 

Jo 

Lt lies in C ([0, T]] L'^{ijj^)) n C ((-oo, Q\x;L'^{ILy x (0, T))) and satisfies the following 
estimates for every t G [0, T] and every x <Q, 

IIC-(-,OL-< IIC.L- + / ll/-v^-I7,,_L-(s)ds, (5.17) 

||C-(a:,.)lk. <^ (\Mu.- + j\\f -^h-U^^_\\^-{s)ds^ . (5.18) 

Once the solutions of (|5.ip - (|5.2p are known it is possible to define the transmission 
conditions on the water- height for the next step (see p.30p ) 

BiX'l+^ := tioC-+i(0,-)) + Jl"+'(0,-). (5.19) 

In the domain x > 0, the situation is slightly different since there are ingoing char- 
acteristics on 7t. So we choose test functions that do not necessarily vanish on the 
boundary and use the transmission condition to prescribe the value of the solution on 
7t. Finally, a solution of (|5.4p satisfies 

Vx e I?(EJ+ X [0,T)), - [C+Adt + t/o • V4x)^+ 

= (C.,x(-,0))^++(Cb,x(0,.))K, + (/-V„.C7;,^+,x)^+, (5.20) 
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where the boundary value (^b is defined on 7^ by 

Cb := —{Bix^-u+}. (5.21) 



Definition 5.5. Let f e L'^{u;t), ?7,, + (= [7^+^) G i2(o^T;V+), G € ^^(tj). 
Assuming that defined by (|5.2ip belongs to L'^{'^t) , we say that C+^^ = C+ G 
L'^{llj^) is a weak solution of Problem (|5.4p z/ (|5.20p holds. 

Using the characteristic method, we have 

Proposition 5.6. Let f, Uf^_^ {— UJ^\^), d and d, be as in Definition 15.51 
There exists a unique weak solution C+^^ G L^{uj^) of (15. 4p . Moreover it is given by 
the characteristic formula: 

C+{x,y,t) = (i{x - UQt,y - vot) + / {f - V h ■ U ^ ^){x - uoS,y - VQS,t - s)ds 

Jo 

if X > uqI, and 

C+{x,y,t) = Cb i y - —x,t - — ] + {f -Vh-U,^^){x-uos,y~vos,t- s)ds, 
\ Uo uoj Jo 

with given by (|5.2ip . if x < uot. 

The solution belongs to C {[0,T]; L'^{io+)) C {[0,+oo)^;L'^{ILy x (0,r))) and sat- 
isfies the following estimates for every t G [0, T] and every x > 0, 

liC+(-,t)L+< \\d\U+uo\\Cbh,+ f \\f-Vh-U„,+ \Uis)ds, (5.22) 

Jo 

\\C+{x,-)K<-(\\dL++uo\\Cb\U,+ f \\f-^h-Uh.+LAs)ds). (5.23) 
"0 V Jo / 

5.3. Well-posedness of the algorithm. First we define a weak formulation 
for the left and right sub-problems at step n of the algorithm. 

Definition 5.7. Let Y = (F,/) e L^i^r) x L^itor), let X, = (C/,,.,,Ci) G 
L^{n) X L^{uj). Forn> 0. 

• Let B1 e L^{0,T;W'). Then X'l+^ = iU^_X-) is a weak solution of Prob- 
lem (|5.ip . (|5.2p i/ it /las regularity (|5.5p - (l5.6p anrf «/ C/^ _ (respectively C,-) is 
a weak solution of ^^(respectively (15. 2p ). 

• Let e L2(-o^7^.)4;/) ^C^r> g l2(7t). T/ien = (C/^_^,C+) is 
a zweafc solution of Problem ()5.3p . ()5.4I1 i/ it /las regularity ()5.5p - (l5.6p anrf i/ 
Uf^ (respectively C+) is a weak solution of ^^(respectively (|5.4p ). 

Then we give a weak formulation for the complete algorithm. 



Definition 5.8. The weak formulation of Algorithm (|3.32p is defined by 



Choose bIx^ £ L^{jt) and G L^{0,T-W) 



Then, for n > 0, 
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• Find weak solution of (|5.ip - (l5.2p and weak solution of 

• Define the transmission conditions for step n + 1 by (j5.12p and (|5.19p . 

Theorem 5.1. FFif/i the hypotheses of Definition \5. i/iere exists a unique weak 
solution XH:'^'^ (respectively Xl+^ ) of Problem dm]) , ([521) (respectively ((0)) . ([5:i)) ). 

Proof. We only prove the result for the left sub-problem, the other one being 
similar. As in the proof of Theorem l4.11 we use a fixed point method. Let us introduce 
the spaces 

£^:=C{[Q,T],H-) L^{Q,T;V-), 

£^:^C{[0,TlL\uj-)) n C((-oo, 0],, (R^ x (0, T)*)) . 

Proposition IS . 21 (respectively Proposition l5.4p defines an afline mapping 
C- 1-^ [//j _ (respectively ^2 : ^ Uf^ _ i-^ (-). 

Setting £^ := £^x£^ — , an application X"^^ is a weak solution of Problem (|5.ip . (|5.2p 
if and only if it is a fixed point in £j^ of the mapping 

T- : (C/,,,_,C-)-(5r(C-),52-(C/;,,_)). 

We now show that T~ has a unique fixed point. Let Xi,X2 € £j^ and let (C/^ _i C-) 
Xi - X2 and C) := T-{Xi) - T-{X2). By linearity (l5l4ll yields: for < t < T, 

f\\^U„,_\\l^is)ds^K f \\U,^J\l^is)ds < ^^{IIC-II^^. + IIC-IIU- 
Jo Jo 

And from Gronwall lemma, we obtain for < t < T, 



\\UH,-\m+ / \\yU,^_rn{s)ds 
Jo 

< ne-^ \ tsnpU-{-,s)\\l + snp\\CA^r)\\l] ■ (5.24) 

[ [0,t] R- J 

Now from ((CT7)l and (|??T5)) . we get 

\\C\\l^{t)+uo\\U^r)K < t\\VU^^_\\l- forO<t<T. (5.25) 

Finally, we endow with the norm ||(C^;i._, C-)llf- •= 



1/2 



sup ||C/^,_(-, s)||^_ + ||VC/^,_||^_ + sup s)\\l + 2«;e«^sup \\Ux, OH^, 

,[0,t] ' [0,t] R- 



With this norm (|5.24p ()5.25p imply that for T' e (0, T] small enough, is contracting 
in £t' ■ This yields the existence of a unique fixed point of in fj^, . Wc obtain the 
result on [0, T] by continuation. □ 

Finally, we can state 

Theorem 5.2. The algorithm (|5.8p is well-defined. 
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Proof. We only have to check that for each step the hypotheses of Theorem lS.ll are 
satisfied. The solutions build at step n have regularity (|5.5p . (|5.6p . We easily 

deduce that defined by belongs to L'^{0,T;W) and B'^+X":+^ defined 

by (|5.19p belongs to L'^{jt)- Thus the hypotheses of Theorem 1 5. II hold for step n + 1. 
□ 

Remark 5.9. Although we do not exhibit a proof here, we are able to establish 
the convergence of the algorithm in some cases. More precisely, if the matrices A 
and B defined by (|3.28p are replaced by diagonal matrices A and B, A being positive 
definite and B being non negative, then the algorithm converges. The proof relies on 
the energy method developed for the Shallow water equations without advection term 
in 125 1. Modifying slightly the proof, we can allow A and B to have non vanishing 
skew- symmetric off-diagonal parts. This generalization still does not cover the situa- 
tion (j3.28p because B has a symmetric non vanishing off- diagonal part. Nevertheless, 
numerical evidences of the convergence of the algorithm are given in the next section. 

6. Numerical results. 

6.1. Numerical scheme in the subdomains. For the numerical applications 
we consider for simplicity a 2 dimensional domain and the related two dimensional 
{x, z) version of the primitive equations p.l3p - ()2.17p . Note that the transmission con- 
ditions (|3.26p - (|3.27p arc independent of the transverse y- variable and are not affected 
by this simplification. 

In this subsection we do not deal with the boundary conditions. Hence the processes 
are the same in both subdomains 0=*= and we restrict ourselves to the subdomain ri+ . 
We first describe the space discretization of the subdomains. We consider a regular 
cartesian grid of nx x nz points and we apply a finite volume method. We introduce 
the horizontal space step Ax and the vertical space step Az. For Euler or Navier- 
Stokes type problems it is well known that a good way to recover some numerical 
stability is to compute velocities and pressure on different cells (see for instance [1] 
and the publications devoted to the so-called C-grids). Here we only deal with the 
horizontal velocity and the water height (depending only depending on x and t) plays 
the role of the pressure. We thus have to introduce two types of finite volume meshes 
- see Figure IHTTl The first one is a 2d finite volume mesh and is related to the compu- 
tation of the velocities. For i = Q...nx — 1 and j = 0...nz we denote I = i -\- jn^. The 
cells of this first mesh wiU be denoted C/ = + (-Aa:;/2, Ax/2) x (-Az/2, Az/2). 
where the points X'j' stand for X'j' = {0,—H) + {iAx,jAz) (they are represented 
by a black circle in Figure 16. ip . The second grid is a Id finite volume mesh de- 
voted to the computation of the water height. The cells of this second mesh will 
be denoted ^+1/2 = Xi^i/2 + (—Ax/2, Ax/2), where the points x.i^i/2 stand for 
Xi+1/2 ~ {i + 1/2) Ax (they are represented by a circle with a number inside in Figure 

EI]). 

Let us now consider the discretization of the equations. Let us start with momentum 
equation (|2.13p . We integrate it on the time-space cell [^^,^^+1] x Cj. We compute 
the interface fluxes at time tf.^1^2 by classical centered formulas. We recover the 
well-known Crank-Nicolson scheme. It is known to be second order accurate and 
conditionally stable in the L°° norm under a CFL type condition on the time step 
Atfe = tk+i — tfc. This strategy is applied for all the velocity nodes such that the 
neighboring nodes are included inside the considered subdomain. The related discrete 
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Fig. 6.1. Space discretization of 



nx-1/2 



(nz+1) nx- 



dx 



nx-1 



relations stand 

At( 1 2 1 2 1 

2 He He e 



1 



-Dxl Cj.k+l 



At r 1 

2 fee 



Re 



1 2 1 



■Vrk 



1 



(6.1) 



At - - 

vi,k+i + -TrS uoD^vi.k+i - —D^vi^k+i - ——D^vi^k+i + -u/,fc+i 
Z [ He He e 

At( 1 2 1 2 1 

vi.k uoD^vi^k - -fr^xVi,k - -f^.DzVi^k + -ui^k 

I I He He e 



(6.2) 



where D^ui^k — {ui+i.k — ui~i,k)/ denotes a classical approximation of the first 
derivative in space in horizontal direction, D^uj^k = («/+!, fc — 2u7_fc + u/_i^fc)/(2Ax) 
and D^uj^k = {'u.i+nx,k — '^uj^k + u/-na;,fc)/(2Ax) denote classical approximations of 
second derivatives in space in horizontal and vertical directions, respectively. 



Let us now consider the mass equation p.l6p . We integrate it on time space cells 
[tk, tk+i] X Ci+i/2 - except for i — where we need to use the transmission conditions. 
We compute the interface fluxes by using explicit upwind formulas. The resulting 
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scheme is known to be first order and also conditionally stable under a CFL type 
condition. The related formula stands 

f At \ At 
AtAz 



where 



—+,n+l 



2 



E 



denotes the discrete mean velocity of the flow along the vertical direction. 

6.2. Numerical discretization near the boundaries. We now have to ex- 
plain how we compute the numerical solution when one of the interfaces of the cell 
belongs to the physical boundaries of the domain or to the fictitious one that is related 
to the domain decomposition method. For all cases we choose to work in the same 
finite volume framework that we use in the interior of the subdomains. 



For the physical boundary conditions (|2.14p we use the ghost cells method. This 
method consists in introducing a fictitious cell along the boundary and then using the 
same finite volume strategy as in the interior domain. For no slip conditions ()2.14p 
we choose the values of the unknowns in the fictitious cell to be equal to their values 
in the neighboring interior cell. 



Let us now focus on the numerical treatment on the cells that are connected with 
the interface F = 9$!+ n - see Fig. 16.11 Here we will use a discrete version of the 
transmission conditions (|3.26p - (|3.27p . This discrete information will be the only data 
that will be transmitted from a subdomain to the other one. Let us first consider the 
mass equation (|2.16p . We integrate it on the cell C3/2 to obtain 



Aa; 



/-+,n 
^-1/ 



At 



/-+,n+l 
^l/2,k 



/•+," + 1 
''-1/2, fc 



—+,n+l 
"0,fc 



0. 



where Ciy'^^^i denotes the water height computed in cell c'^^^ at time and for 
iteration n + 1 of the algorithm. The quantities C-i/t]: ^^'^ ^o^fc"^^ have to be consid- 
ered as unknown quantities since the corresponding cells are not included in f2+. We 
will use the transmission conditions p.32p to evaluate them. Hence we obtain thanks 

to (p:^ 



''l/2,fc+l ^ I ^ 



uoAt\ 
Ax I 



^l/2,fe 



At 



Ax 



h,k 



where has been computed in f7 
given by 



during the previous Schwarz iteration and is 



— ,n 
k 
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The basic idea is the same for the momentum equation (|2.13p . Here we integrate the 
equation on the semi-ceh (7/ = + (0, Ax/2) x (— Az/2, Az/2). for / ~ jn^ with 
j = 0, ...,nz. We obtain 



AxAz -^n+l +.n+lN 

{/ N I +.n+l/ N \ 
YI 2 "/,fe+l/2(0 I 

-]^(^ ^ a^-Wi/2(0j 

^ I ^l/2,k ^ '^l/2,fc+l >+.n+l \ I 

"^i^r2 I 2 ''-i/2,fc+i/2 I r 

2 Re' 2 2 £ 2 



where Z (respectively r) denotes quantities that are evaluated on the left (respec- 
tively right) boundary of the cell Cf . Hence quantities ^^'^+1/2^^)^ ^^""z 'fe+i/2(') ^"^"^ 
''-i/2^fc+i/2 ^"^^ unknown quantities since they involve quantities that are computed 
outside the domain 17+ . Here also we use transmission conditions (I3.32p and we obtain 
thanks to ([SlSOl) 

, 1/ ^M:'ir) 1 2D^uZ':,\r) 1 ^_„^, 1 



"0 1 1 — ^^W/.fc+l - 



At 21 Ax i?e Ax Re' ^ ■''*=+^ e 

+,ri+l . +,n+l _ I ^ a\ - + ."+! _ ;^- + .ri+l 




1 f 1 2jJ..;:^'+^(r) 1 ^ 1 ^ 

""^^ -Re A^ i?^^^"^^^ -T'^^^ 

I,k -L ^1/2, fc ^l/2,fe+l 



At Fr^ Ax 



(6.3) 



where j j. has been computed in fl during the previous Schwarz iteration and is 
deduced from relation p.3ip 



— ,n . — ,n — ,n I — ,n 



Note that i3"'" , is known since it has been computed in fl^ at iteration n — 1 and has 
been transmitted to the domain before iteration n. Same type of computations 
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for the transverse component of the velocity lead to the following scheme 







At 






Ax 


Re 


Ax 


+ ,n+l 
'I,k+1 


^ +.n+l 


1 

2 



Re' '^^^ £ r 



1 



A2;^+J.fc ' Ax ] ^^^^^ 




2r°^^ — i?; — — ite^^^"^.^ + T'^^ 

where S^'" has been computed in during the previous Schwarz iteration and is 
deduced from relation (|3.3ip 

— ,n . — .n — I — 

+j,fc ~ y/e 2 y/e 2 

The derivation of the discrete boundary condition in is based on the same type 
of computations. Note that only the components of the velocity are concerned by the 
transmission problem in fi". 

6.3. Numerical optimization of the transmission conditions. In this sec- 
tion we are interested in the optimization the transmission conditions p.30p with 
respect to the free parameters a and f3. To optimize the conditions means that we 
choose parameters a and /3 such that the Schwarz waveform relaxation algorithm 
p.32p reaches a given error for as small as possible number of iterations. The analyti- 
cal solution of this problem is quite complex in the considered framework and we only 
present here a numerical strategy to reach the optimum. In the simpler case of a ID 
advection diffusion equation a complete solution of the related optimization problem 
is given in [llj . 

We consider a test case for which all the initial data (velocities and perturbation 
of the water height) are taken equal to zero. We initialize the algorithm (|3.1[) with 
random boundary conditions on the interface and we study the convergence of the so- 
lution towards the analytical ones. This test is quite classical to study the convergence 
of a domain decomposition algorithm. It is interesting since the initial quantities do 
contain all frequencies. In all the computations the physical parameters Re and Fr 
are taken equal to one but the Rossby number e remains free. For a given value of e 
we apply the transmission conditions (|3.30p for several values of the parameters a and 
f3 and we compare the error between the computed and the analytical solutions 
after a given number of iterations. It allows us to find an optimal pair {uopt, Popt) 
that minimizes this error. This first study exhibit that the influence of the parameter 
P is quite small. In the following this parameter will be kept equal to its theoretical 
value (j3.29p . In a second step we study the dependency of the optimal parameter 
Oiopt with respect to the Rossby number e. The results are presented in Fig. 16.21 for 
different values of e. We found that this optimized parameter does depend on e in a 
nontrivial way. 

We now present the evolution of the error on the computed solution as a function of 
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1 I ■ — ' — ■ — ■ — ■ — ' — ' — ■ — ' — 

0.0001 0.001 0.01 0.1 

Fig. 6.2. Quotient Qopt/'^Tay between the numerically optimized parameter and the Taylor 
approximation parameter as a function of the Rossby number e (in Log scale) 



the number of iterations of the Schwarz waveform relaxation algorithm (|3.ip in both 
cases a = Uopt and a = axay ■ in Fig. 16.31 we present the results for two different 
values of the Rossby number e : e ~ 10^^ and e ~ 10^^. The curves (Log of the 
error) all look like straight lines, at least after a sufficiently large number of iterations. 
The method appears to be more efficient when the Rossby number is smaller since the 
error decreases much faster in the case e = 10^^ - Fig. 16.31 on the left. This result is 
consistent with the previous theoretical study that is based on an asymptotic analysis 
in e. We also observe that for a given value of e the curves look similar for both 
optimized and Taylor approximation parameters even if the error decreases faster for 
the optimal value aopt- Moreover let us observe that to reach an error of 10"'* (that 
is enough for the applicability of the Schwarz waveform relaxation algorithm) both 
algorithms (with optimized or Taylor approximation parameter) need a very close 
number of iterations. 




5 ID 15 20 5 ID 15 20 



Fig. 6.3. Log of the error on the computed solution as a function of the number of Schwarz 
iterations for Rossby number e = 10"'^ (l^ft) and e = 10^^ (right) and for a random initial guess 
using the Taylor approximation parameter QTay (^p) o-iT-d the optimized one ctopt (down) 
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We compute the same test with Rossby number e = 10^^ but with a sinusoidal ini- 
tial guess (instead of the random ones) for the transmission conditions. We consider 
two different sinusoids with one or ten periods in the space-time considered interval 
and we use Taylor approximation parameters axay and Pray In Fig- 16.41 the results 
appears to be much better for the low frequency sinusoid as for high frequency one. 
The results for the high frequency sinusoid look similar to the results that were ob- 
tained with the random initial guess. It follows that the method is particularly well 
adapted to low frequency signals : the relative error is smaller than 10^^ after only 
two iterations. 




5 10 15 20 



Fig. 6.4. Log of the error on the computed solution as a function of the number of Schwarz 
iterations for Rossby number e = 10~^ and with optimal parameter for a low frequency signal ( down) 
and for a high frequency signal (middle) and for a random signal (up) 

6.4. Numerical application. In this section we consider the case of a flow with 
a constant positive background velocity uq = l.m/s and an initial local decreasing 
step on the water height. We choose the Rossby number e equal to 10"'^. We choose 
nx = 40, nz — 10 and nt = 40 in order to ensure the CFL condition. We present the 
initial solution and the solution computed at final time T — 1.3s after 20 iterations by 
the proposed Schwarz waveform relaxation algorithm in Fig. 16.61 The 2d horizontal 
velocity vector field (u, v) is presented in the 2d vertical domain (in the {x, z) plane) 
which is occupied by the flow. A horizontal vector denotes a velocity which is coUinear 
to the x-direction and a vertical one denotes a velocity which is coUinear to the y- 
direction. Since we consider the linearized version of the equations the step just moves 
without deformation from the left to the right of the domain. Since the Coriolis effect 
is dominant we observe the formation of a transverse jet which moves with the step. 
Another consequence of the Coriolis effect is the formation of a stationary eddy at the 
initial location of the step. We now compare the solution that is computed on the 
whole domain with the solution that is obtained by considering the presented domain 
decomposition strategy. In Fig. 16.71 we present the evolution of the relative error 
between the two solutions versus the number of considered iterations. It exhibits the 
fast convergence of the algorithm for such a case. After two iterations the relative 
error is around 10~^ and it reaches the factor 10^^" after eight iterations. 

7. Conclusion. We presented in this article a new domain decomposition method 
for the viscous primitive equations. It involves a Schwarz waveform relaxation type al- 
gorithm with approximated transmission conditions for which we proved well-posedness 
We presented a numerical optimization of the transmission conditions and we study 
the speed of convergence of the algorithm for several test cases. Academic numerical 
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Fig. 6.5. Water height and velocity field at initial time 




Fig. 6.6. Water height and velocity field at final time 



applications were presented. In forthcoming papers we plan to prove the convergence 
of the algorithm and we want to present oceanographic configurations and to increase 
the efficiency of the algorithm by deriving more complex transmission conditions based 
on another asymptotic regime that corresponds to quasi-geostrophic flows. 
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Fig. 6.7. Log of the relative error on the solution computed by using the Schwarz waveform 
relaxation algorithm versus number of iterations 
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